Código
## primeras 6 líneas
head(d) y x1 x2 x3 x4 x5
1 2625 20 50 60 B 0
2 2928 20 75 40 A 0
3 2823 15 50 60 B 1
4 2824 15 75 60 C 0
5 2604 20 60 60 B 0
6 2770 30 60 80 C 1
Sesión #3: Ajuste y Validación de un Modelo de RLM
En esta práctica ajustaremos y validaremos un modelo de RLM, y finalmente realizaremos predicción.
Pollos Riko Riko (PRR) es la empresa líder en productos avícolas de la Región Caribe. Su centro de operaciones, ubicado en Sabanagrande, Atlántico, produce 4 tipos de producto: (i) pollo entero, (ii) bandejas de pechuga entera, (iii) bandejas de pernil y (iv) bandejas de alas. El precio promedio de venta de cada producto, por libra, es $8000, $4300, $3400 y $2900, respectivamente. Se sabe que la participación de cada producto en las ventas totales de la compañía es \(0<p_j<1\) conocido, \(j=1,2,3,4\). Por supuesto, \(p_1+p_2+p_3+p_4 = 1\).
Los animales se sacrifican luego de 40 días de ser alimentados con una dieta balanceada que incluye nutrientes especiales (variable \(x_1\) en gramos/día), agua (variable \(x_2\) en ml/día) y forraje (variable \(x_3\) en gramos/día), además de la raza (variable \(x_4\) con niveles A, B y C) y el hecho de que sean expuestos a una luz especial durante la noche (variable \(x_5\) con niveles 0: No y 1: Si). Actualmente, el peso promedio de un pollo que crece en las instalaciones de la compañía está en el intervalo (2400, 2800) gramos, con una confianza del 95%.
Con miras a aumentar la eficiencia de la planta,1 PRR ha decidido aumentar el peso de los animales antes de su sacrificio. Para ello, decide realizar un experimento en el que a 100 grupos de 100 animales (i.e., lote) se les proporciona la dieta balanceada y se cuantifica, al final del tiempo de engorde, el peso promedio alcanzado (variable respuesta \(Y\)).
Para leer los datos hacemos:
Las primeras 6 filas de la base de datos d son
## primeras 6 líneas
head(d) y x1 x2 x3 x4 x5
1 2625 20 50 60 B 0
2 2928 20 75 40 A 0
3 2823 15 50 60 B 1
4 2824 15 75 60 C 0
5 2604 20 60 60 B 0
6 2770 30 60 80 C 1
Analicemos inicialmente la correlación entre las variables disponibles:
## matriz de correlación
par(mfrow = c(1,1), mar = c(.1, .1, .1, .1))
qgraph(cor(d[, -c(5, 6)]), graph = "cor", layout = "spring",
sampleSize = nrow(d),
legend.cex = 1, alpha = 0.05)Numéricamente, la matriz de correlación es
## matriz de correlación
cor(d[, -c(5, 6)]) y x1 x2 x3
y 1.0000000 0.18310750 0.41359388 -0.4468526
x1 0.1831075 1.00000000 0.05908606 -0.0254010
x2 0.4135939 0.05908606 1.00000000 -0.1144729
x3 -0.4468526 -0.02540100 -0.11447293 1.0000000
Estos resultados indican que la correlación entre \(y\) y \(x_1\) es 0.183, entre \(y\) y \(x_2\) es 0.413, y entre \(y\) y \(x_3\) es -0.447. En cuanto a la correlaciones entre las variables independientes, podemos concluir que estas son pequeñas, lo cual sugiere que, efectivamente, \(x_1, x_2\) y \(x_3\) son independientes.
También podemos representar las correlaciones y la distribución de cada variable en un gráfico de dispersión/correlación:
## gráfico de dispersión/correlación
GGally::ggpairs(d) + theme_minimal()Más información, aquí.
El gráfico 3D entre \(x_1\), \(x_2\) y \(y\) sería:
fig <- plot_ly(d,
x = ~x1,
y = ~x2,
z = ~y,
text = ~rownames(d),
color = '#BF382A')
fig <- fig %>% add_markers()
fig <- fig %>% layout(title = '\n y vs. (x1, x2)',
scene = list(xaxis = list(title = 'x1'),
yaxis = list(title = 'x2'),
zaxis = list(title = 'y')))
figEl modelo ajustado es:
## full MLR model
fit <- lm(y ~ ., data = d)
summary(fit)
Call:
lm(formula = y ~ ., data = d)
Residuals:
Min 1Q Median 3Q Max
-222.022 -51.849 7.082 57.180 187.069
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2630.7862 68.4619 38.427 < 2e-16 ***
x1 2.0990 1.2563 1.671 0.098135 .
x2 4.1656 0.8549 4.873 4.51e-06 ***
x3 -2.5947 0.5130 -5.058 2.12e-06 ***
x4B -86.5364 21.6830 -3.991 0.000131 ***
x4C -10.0919 22.2276 -0.454 0.650868
x51 74.8475 17.8246 4.199 6.13e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 84.65 on 93 degrees of freedom
Multiple R-squared: 0.5742, Adjusted R-squared: 0.5468
F-statistic: 20.9 on 6 and 93 DF, p-value: 2.237e-15
Los intervalos de confianza del 95% para los coeficientes pueden obtenerse a través de la función confint.default() haciendo
## 95% CI para los coeficientes
confint.default(fit) 2.5 % 97.5 %
(Intercept) 2496.6034309 2764.969042
x1 -0.3633595 4.561367
x2 2.4901124 5.841087
x3 -3.6001142 -1.589352
x4B -129.0344167 -44.038429
x4C -53.6572601 33.473418
x51 39.9118304 109.783094
En términos generales, el concepto de multicolinealidad es sinómino de redundancia en las variables independientes.
Uno de los supuestos fuertes del modelo de RLM es que las variables \(X_1, X_2,\ldots,X_n\) son independientes. Cuando esto no ocurre, los estimadores de \({\beta} = (\beta_1,\beta_2,\ldots,\beta_k)\) tienen propiedades distintas a las ya conocidas.
Desde el punto de vista formal, la existencia de multicolinealidad puede probarse utilizando el ill-condition number (ICN)
## ICN
kappa(fit) [1] 342.3963
Teniendo en cuenta que el ICN es \(>30\), aparentemente, existe multicolinealidad entre \(x_1, x_2\) y \(x_3\).
Ahora, si estamos interesados en determinar cuál de la(s) variable(s) independiente(s) con mayor grado de colinealidad, utilizamos el VIF:
## VIF
car::vif(fit) GVIF Df GVIF^(1/(2*Df))
x1 1.082171 1 1.040275
x2 1.047698 1 1.023571
x3 1.028809 1 1.014302
x4 1.082995 2 1.020132
x5 1.108078 1 1.052653
Al analizar el VIF, es claro que ningún valor es \(>5\).
Recientemente, se han implementado otras pruebas de multicolinealidad en el paquete mctest:
## otras pruebas de multicolinealidad
mctest(fit)$odiags results detection
Determinant 0.5730669 0
Farrar Chi-Square 53.5410606 1
Red Indicator 0.1786286 0
sum of Lambda Invers 7.3543905 0
Theil Indicator -1.9206744 0
Condition Number 22.2467502 0
De acuerdo con estos resultados, podría existir multicolinealidad. En particular, aparece el valor de 1 en la segunda entrada de la columna detection.
La validación de los supuestos puede hacerse via valores \(p\) como se muestra a continuación.
Normalidad
## prueba de Normalidad
r <- rstudent(fit)
shapiro.test(r)$p.value[1] 0.7810988
Como el valor \(p\) es \(>0.05\), entonces los errores del modelo siguen una distribución Normal.
Varianza constante
## prueba de varianza constante
car:::ncvTest(fit)$p[1] 0.7729436
Como el valor \(p\) es \(>0.05\), podemos concluir que los errores tienen varianza constante.
Independencia
## prueba de independencia
car:::durbinWatsonTest(fit)$p[1] 0.436
Este resultado indica que los errores del modelo ajustado son independientes.
Para identificar outliers usamos los residuales estudentizados del modelo:
## outliers?
res <- which(r > 3 | r < -3)
ifelse(length(res) == 0, 0, length(res))[1] 0
Otra forma de detectar outliers es a través de la prueba de Bonferroni. Esta prueba está implementada en la función outlierTest del paquete car. En nuestro ejemplo, procedemos de la siguiente manera:
## prueba Bonferroni para outliers
outlierTest(fit, n.max = sqrt(NROW(d)))No Studentized residuals with Bonferroni p < 0.05
Largest |rstudent|:
rstudent unadjusted p-value Bonferroni p
31 -2.840257 0.0055492 0.55492
Este resultado indica que hay 0 outliers en los datos.
Para identificar este tipo de observaciones, utilizamos la distancia de Cook.
En R procedemos de la siguiente manera:
## gráfico de la distancia de Cook
plot(fit, which = 4, las = 1)Este resultado indica que las observaciones 20, 30 y 31 podrían considerarse influenciales.
Conclusión: Los supuestos de Normalidad, independencia y varianza constante de los errores se cumplen. Por lo tanto, el modelo es válido para predecir. Además, parecen no existir outliers en los datos.
Si fuese de interés determinar el peso promedio de los pollos con en las condiciones
\[ \mathbf{x}_0= (28, 65, 70, \text{A}, 1) \]
procedemos de la siguiente forma:
## x0
x0 <- data.frame(x1 = 28, x2 = 65, x3 = 70, x4 = 'A', x5 = "1")
x0 x1 x2 x3 x4 x5
1 28 65 70 A 1
## estimación
predict(fit, newdat = x0) 1
2853.538
Por lo tanto, \(\widehat{E[Y|\mathbf{x}_0]} = 2778.7\).
interval a predict():El intervalo de confianza del 95% es
## confidence interval
predict(fit, newdat = x0, interval = 'confidence') fit lwr upr
1 2853.538 2809.43 2897.647
## prediction interval
predict(fit, newdat = x0, interval = 'prediction') fit lwr upr
1 2853.538 2679.753 3027.324
Así, el peso del próximo pollo engordado en las condicione \(\mathbf{x}_0= (28, 65, 70, \text{A}, 1)\) será
\[ Y | \mathbf{x}_0 \in (2679.7, 3027.3). \]
Para más información sobre el cálculo de los intervalos de confianza y predicción a partir de un modelo de RLM ajustado, se recomienda consultar ?predict.lm.
Fecha de entrega: Abril 12, 2024.
En kilogramos:
ancho <- 14
largo <- 140
area <- ancho*largo
n_galpones <- 200
total_n_pollos = area*n_galpones*8
print(total_n_pollos*mean(d$y)/1000)[1] 8689511
NO podemos decir que NO tiene sentido, de hecho utilizar la luz es una varaible significativa según el modelo construido, además es la tercera varaible más importante y la utilizar la luz especial supone un aumento de 74.84 g en el peso esperado de los pollos luego de la etapa engorde, en comparación de no utilizar la luz.
La variable con mayor importancia es x3, la cantidad de forraje que le dan a los pollos en gramos/días. Si bien su magnitud es de -2.5, es la variable con menos error estándar, lo que añade relevancia a su resultado.
B y determine el peso promedio esperado cuando \(\mathbf{x}_0= (28, 65, 70, \text{B}, 1)\).Recomendaría la raza A, ya que existe sufiente evidencia estadística para afirmar que esta raza representa un aumento de 86 g del valor esperado en comparación con la raza B. Y aunque no se puede decir lo mismo en compración con la raza C, la media de la raza A en los datos de muestra es 10 g superior, pero de no poder utilizar la raza A, los resultados sugieren que no hay una diferencia estadística con la raza C.
mean(d[d$x4=="A",]$y)[1] 2822.815
mean(d[d$x4=="C",]$y)[1] 2810.375
No es posible hablar de uniformidad en el peso independiente de la raza, como se menciona anteriormente, por lo menos la raza B tiene una diferencia estadística con las otras razas.
Para las condiciones dadas, se espera un peso promedio estimado para un galpón de 2764 gramos:
x0 <- data.frame(x1 = 28, x2 = 65, x3 = 70, x4 = 'B', x5 = "1")
predict(fit, newdat = x0) 1
2767.002
x1 <- data.frame(x1 = 28, x2 = 65, x3 = 70, x4 = 'B', x5 = "1")
peso_esperado <- predict(fit, newdat = x1)[[1]]
total_n_pollos_1glp <- area*8
# Total número de pollos en 40 días
peso_total_lb_1glp <- total_n_pollos_1glp*peso_esperado/453
producto1 <- 0.1*8000*peso_total_lb_1glp
producto2 <- 0.3*4300*peso_total_lb_1glp*0.40
producto3 <- 0.45*3400*peso_total_lb_1glp*0.30
producto4 <- 0.15*2900*peso_total_lb_1glp*0.15
#total ingresos mensuales para estos productos
ingresos <- (producto1+producto2+producto3+producto4)*30/40Con esto, por galpón se espera una utilidad de 122 millones de pesos, sin embargo, esto es sólo un escenario y existen otras configuraciones que podrían generar más utilidad. Por ejemplo, cambiar la raza de los pollos utilizados.
Para más información sobre el modelo de RLM se sugiere consultar el Capítulo 3 del texto Modelos de Regresión: Una aproximación práctica con R.
Esto se refiere a que, al final del período de engorde, el animal pese más de lo que pesa en las condiciones actuales.↩︎